import numpy as np

def newton_generation(initial_value, e, M, tolerance):
    """
    牛顿迭代求解开普勒方程 E - e sin E = M
    :param initial_value: 开始迭代的值，可以为任意值，下面代码中取了上一个时刻的E为初始值
    :param e: 轨道根数
    :param M: 同上
    :param tolerance: 迭代到误差小于什么时候停止，下面代码中取了1e-6
    :return: 计算结果
    """
    difference = 2 * tolerance
    old_value = initial_value
    new_value = initial_value
    while difference > tolerance:
        new_value = old_value - (old_value - e * np.sin(old_value) - M) / (1 - e * np.cos(old_value))
        difference = abs(old_value - new_value)
        old_value = new_value

    return new_value

def rk4(x, t, h, f):
    """
    rk4数值法求解微分方程 dx/dt = f(x,t)
    :param x:
    :param t:
    :param h:
    :param f:
    :return:
    """
    k1 = h * f(x, t)
    k2 = h * f(x + 0.5 * k1, t + 0.5 * h)
    k3 = h * f(x + 0.5 * k2, t + 0.5 * h)
    k4 = h * f(x + k3, t + h)
    return x + (k1 + 2 * k2 + 2 * k3 + k4) / 6